Decorrelation-Based Feature Discovery: COVID_19

1 Effect of UPSTM-Based Decorrelation on Feature Discovery

1.0.1 Loading the libraries

library("FRESA.CAD")
library(readxl)
library(igraph)
library(umap)
library(tsne)
library(entropy)
library(TH.data)
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)

1.1 Material and Methods

About Dataset Coronaviruses are a broad family of viruses that have been linked to illnesses ranging from the common cold to more serious conditions such as Middle East Respiratory Syndrome (MERS) and Severe Acute Respiratory Syndrome (SARS) (SARS. In 2019, a new coronavirus (COVID-19) was discovered in Wuhan, China.

Sometimes, When PCR test resources are scarce and antigen test kits are inaccurate, clinicians look for alternate COVID-19 testing methods that can be completed in a day and handle thousands of samples. COVID-19 virus proteins should be lacking in normal people’s saliva. Some distinct proteins may be produced in response to COVID-19 infection and can be utilized as a signature to identify potentially infected people. Therefore, protein profiles in a patient’s saliva can indicate that he or she is infected with COVID-19.

Mass spectrometry is a method for determining the protein composition of a material. Saliva samples from hundreds of patients were studied in this dataset. So, a machine learning specialist is approached and charged with developing a machine learning model that can identify who COVID-19 infected since the PCR test cannot interpret the result completely.

1.2 The Data

https://www.kaggle.com/datasets/kerneler/saliva-testing-dataset?select=COVID-19_MS_dataset_train.csv

COVID_19_MS <- read.csv("~/GitHub/LatentBiomarkers/Data/COVID19_MS/COVID-19_MS_dataset_train.csv")

colnames(COVID_19_MS) <- str_replace_all(colnames(COVID_19_MS),"\\.","_")
COVID_19_MS$Person_ID <- NULL
sampleID <- unique(COVID_19_MS$Sample_ID)
spectraID <- colnames(COVID_19_MS)[str_detect(colnames(COVID_19_MS),"X")]

avgCOVID19 <- NULL
class <- COVID_19_MS[!duplicated(COVID_19_MS$Sample_ID),"PCR_result"]
for (id in sampleID)
{
  avgCOVID19 <- rbind(avgCOVID19,apply(COVID_19_MS[COVID_19_MS$Sample_ID %in% id,spectraID],2,mean))
}

avgCOVID19 <- as.data.frame(avgCOVID19)
rownames(avgCOVID19) <- sampleID
avgCOVID19$class <- 1*(class=="pos")

pander::pander(table(avgCOVID19$class))

0	1
37	109

1.2.0.1 Standarize the names for the reporting

studyName <- "COVID19_MS"
dataframe <- avgCOVID19
outcome <- "class"
thro <- 0.80
TopVariables <- 10
cexheat = 0.25

1.3 Generaring the report

1.3.1 Libraries

Some libraries

library(psych)
library(whitening)
library("vioplot")
library("rpart")

1.3.2 Data specs

pander::pander(c(rows=nrow(dataframe),col=ncol(dataframe)-1))

rows	col
146	2715

pander::pander(table(dataframe[,outcome]))

0	1
37	109

varlist <- colnames(dataframe)
varlist <- varlist[varlist != outcome]

largeSet <- length(varlist) > 1500 

1.3.3 Scaling the data

Scaling and removing near zero variance columns and highly co-linear(r>0.99999) columns

  ### Some global cleaning
  sdiszero <- apply(dataframe,2,sd) > 1.0e-16
  dataframe <- dataframe[,sdiszero]

  varlist <- colnames(dataframe)[colnames(dataframe) != outcome]
  tokeep <- c(as.character(correlated_Remove(dataframe,varlist,thr=0.99999)),outcome)
  dataframe <- dataframe[,tokeep]

  varlist <- colnames(dataframe)
  varlist <- varlist[varlist != outcome]
  
  iscontinous <- sapply(apply(dataframe,2,unique),length) >= 5 ## Only variables with enough samples



dataframeScaled <- FRESAScale(dataframe,method="OrderLogit")$scaledData

1.4 The heatmap of the data

numsub <- nrow(dataframe)
if (numsub > 1000) numsub <- 1000


if (!largeSet)
{

  hm <- heatMaps(data=dataframeScaled[1:numsub,],
                 Outcome=outcome,
                 Scale=TRUE,
                 hCluster = "row",
                 xlab="Feature",
                 ylab="Sample",
                 srtCol=45,
                 srtRow=45,
                 cexCol=cexheat,
                 cexRow=cexheat
                 )
  par(op)
}

1.4.0.1 Correlation Matrix of the Data

The heat map of the data

if (!largeSet)
{

  par(cex=0.6,cex.main=0.85,cex.axis=0.7)
  #cormat <- Rfast::cora(as.matrix(dataframe[,varlist]),large=TRUE)
  cormat <- cor(dataframe[,varlist],method="pearson")
  cormat[is.na(cormat)] <- 0
  gplots::heatmap.2(abs(cormat),
                    trace = "none",
  #                  scale = "row",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "Original Correlation",
                    cexRow = cexheat,
                    cexCol = cexheat,
                     srtCol=45,
                     srtRow=45,
                    key.title=NA,
                    key.xlab="|Pearson Correlation|",
                    xlab="Feature", ylab="Feature")
  diag(cormat) <- 0
  print(max(abs(cormat)))
}

1.5 The decorrelation

DEdataframe <- IDeA(dataframe,verbose=TRUE,thr=thro)
#> 
#>  X3404_08 X3325_02 X5072_49 X348_08 X4279_48 X11050_18 
#>   X244_38   X245_48   X246_44   X246_69   X247_36   X247_99 
#> 0.6583199 0.6587237 0.4390145 0.3945880 0.6817447 0.4886914 
#> 
#>  Included: 2476 , Uni p: 6.058158e-05 , Base Size: 43 , Rcrit: 0.312781 
#> 
#> 
 1 <R=0.994,thr=0.950>, Top: 12< 1 >[Fa= 12 ]( 12 , 22 , 0 ),<|><>Tot Used: 34 , Added: 22 , Zero Std: 0 , Max Cor: 0.975
#> 
 2 <R=0.975,thr=0.950>, Top: 2< 3 >[Fa= 14 ]( 2 , 5 , 12 ),<|><>Tot Used: 41 , Added: 5 , Zero Std: 0 , Max Cor: 0.968
#> 
 3 <R=0.968,thr=0.950>, Top: 1< 1 >[Fa= 15 ]( 1 , 1 , 14 ),<|><>Tot Used: 43 , Added: 1 , Zero Std: 0 , Max Cor: 0.949
#> 
 4 <R=0.949,thr=0.900>, Top: 26< 2 >[Fa= 32 ]( 26 , 42 , 15 ),<|><>Tot Used: 98 , Added: 42 , Zero Std: 0 , Max Cor: 0.945
#> 
 5 <R=0.945,thr=0.900>, Top: 3< 2 >[Fa= 35 ]( 3 , 4 , 32 ),<|><>Tot Used: 105 , Added: 4 , Zero Std: 0 , Max Cor: 0.900
#> 
 6 <R=0.900,thr=0.800>, Top: 74< 1 >[Fa= 84 ]( 73 , 141 , 35 ),<|><>Tot Used: 287 , Added: 141 , Zero Std: 0 , Max Cor: 0.856
#> 
 7 <R=0.856,thr=0.800>, Top: 8< 1 >[Fa= 90 ]( 8 , 8 , 84 ),<|><>Tot Used: 299 , Added: 8 , Zero Std: 0 , Max Cor: 0.800
#> 
 8 <R=0.800,thr=0.800>
#> 
 [ 8 ], 0.7956086 Decor Dimension: 299 Nused: 299 . Cor to Base: 203 , ABase: 2476 , Outcome Base: 0 
#> 
varlistc <- colnames(DEdataframe)[colnames(DEdataframe) != outcome]

pander::pander(sum(apply(dataframe[,varlist],2,var)))

23332

pander::pander(sum(apply(DEdataframe[,varlistc],2,var)))

4516

pander::pander(entropy(discretize(unlist(dataframe[,varlist]), 256)))

0.132

pander::pander(entropy(discretize(unlist(DEdataframe[,varlistc]), 256)))

0.449

varratio <- attr(DEdataframe,"VarRatio")

pander::pander(tail(varratio))

La_X3405_23	La_X3396_15	La_X570_46	La_X3283_6	La_X547_98	La_X3333_73
0.0307	0.0273	0.024	0.0225	0.0212	0.0114

1.5.1 The decorrelation matrix

if (!largeSet)
{

  par(cex=0.6,cex.main=0.85,cex.axis=0.7)
  
  UPLTM <- attr(DEdataframe,"UPLTM")
  
  gplots::heatmap.2(1.0*(abs(UPLTM)>0),
                    trace = "none",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "Decorrelation matrix",
                    cexRow = cexheat,
                    cexCol = cexheat,
                   srtCol=45,
                   srtRow=45,
                    key.title=NA,
                    key.xlab="|Beta|>0",
                    xlab="Output Feature", ylab="Input Feature")
  
  par(op)
  
  
  
}

1.5.2 Formulas Network

Displaying the features associations

par(op)


#if ((ncol(dataframe) < 1000) && (ncol(dataframe) > 10))
#{

#  DEdataframeB <- ILAA(dataframe,verbose=TRUE,thr=thro,bootstrap=30)

  transform <- attr(DEdataframe,"UPLTM") != 0
  tnames <- colnames(transform)
  colnames(transform) <- str_remove_all(colnames(transform),"La_")
  transform <- abs(transform*cor(dataframe[,rownames(transform)])) # The weights are proportional to the observed correlation
  
  
  VertexSize <- attr(DEdataframe,"fscore") # The size depends on the variable independence relevance (fscore)
  names(VertexSize) <- str_remove_all(names(VertexSize),"La_")
  VertexSize <- 10*(VertexSize-min(VertexSize))/(max(VertexSize)-min(VertexSize)) # Normalization

  VertexSize <- VertexSize[rownames(transform)]
  rsum <- apply(1*(transform !=0),1,sum) + 0.01*VertexSize + 0.001*varratio[tnames]
  csum <- apply(1*(transform !=0),2,sum) + 0.01*VertexSize + 0.001*varratio[tnames]
  
  ntop <- min(10,length(rsum))


  topfeatures <- unique(c(names(rsum[order(-rsum)])[1:ntop],names(csum[order(-csum)])[1:ntop]))
  rtrans <- transform[topfeatures,]
  csum <- (apply(1*(rtrans !=0),2,sum) > 1)
  rtrans <- rtrans[,csum]
  topfeatures <- unique(c(topfeatures,colnames(rtrans)))
  print(ncol(transform))
#> [1] 299
  transform <- transform[topfeatures,topfeatures]
  print(ncol(transform))
#> [1] 20
  if (ncol(transform)>100)
  {
    csum <- (apply(1*(transform !=0),2,sum) > 1) & (apply(1*(transform !=0),1,sum) > 1)
    transform <- transform[csum,csum]
    print(ncol(transform))
  }

    if (ncol(transform) < 150)
    {

      gplots::heatmap.2(transform,
                        trace = "none",
                        mar = c(5,5),
                        col=rev(heat.colors(5)),
                        main = "Red Decorrelation matrix",
                        cexRow = cexheat,
                        cexCol = cexheat,
                       srtCol=45,
                       srtRow=45,
                        key.title=NA,
                        key.xlab="|Beta|>0",
                        xlab="Output Feature", ylab="Input Feature")
  
      par(op)
      VertexSize <- VertexSize[colnames(transform)]
      gr <- graph_from_adjacency_matrix(transform,mode = "directed",diag = FALSE,weighted=TRUE)
      gr$layout <- layout_with_fr
      
      fc <- cluster_optimal(gr)
      plot(fc, gr,
           edge.width = 2*E(gr)$weight,
           vertex.size=VertexSize,
           edge.arrow.size=0.5,
           edge.arrow.width=0.5,
           vertex.label.cex=(0.15+0.05*VertexSize),
           vertex.label.dist=0.5 + 0.05*VertexSize,
           main="Top Feature Association")
      
    }

par(op)

1.6 The heatmap of the decorrelated data

if (!largeSet)
{

  hm <- heatMaps(data=DEdataframe[1:numsub,],
                 Outcome=outcome,
                 Scale=TRUE,
                 hCluster = "row",
                 cexRow = cexheat,
                 cexCol = cexheat,
                 srtCol=45,
                 srtRow=45,
                 xlab="Feature",
                 ylab="Sample")
  par(op)
}

1.7 The correlation matrix after decorrelation

if (!largeSet)
{

  cormat <- cor(DEdataframe[,varlistc],method="pearson")
  cormat[is.na(cormat)] <- 0
  
  gplots::heatmap.2(abs(cormat),
                    trace = "none",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "Correlation after ILAA",
                    cexRow = cexheat,
                    cexCol = cexheat,
                     srtCol=45,
                     srtRow=45,
                    key.title=NA,
                    key.xlab="|Pearson Correlation|",
                    xlab="Feature", ylab="Feature")
  
  par(op)
  diag(cormat) <- 0
  print(max(abs(cormat)))
}

1.8 U-MAP Visualization of features

1.8.1 The UMAP on Raw Data

  classes <- unique(dataframe[1:numsub,outcome])
  raincolors <- rainbow(length(classes))
  names(raincolors) <- classes
  topvars <- univariate_BinEnsemble(dataframe,outcome)
  lso <- LASSO_MIN(formula(paste(outcome,"~.")),dataframe,family="binomial")
  topvars <- unique(c(names(topvars),lso$selectedfeatures))
  pander::pander(head(topvars))

X301_02, X520_38, X3283_6, X3354_77, X3439_22 and X3540_88

#  names(topvars)
#if (nrow(dataframe) < 1000)
#{
  datasetframe.umap = umap(scale(dataframe[1:numsub,topvars]),n_components=2)
#  datasetframe.umap = umap(dataframe[1:numsub,varlist],n_components=2)
  plot(datasetframe.umap$layout,xlab="U1",ylab="U2",main="UMAP: Original",t='n')
  text(datasetframe.umap$layout,labels=dataframe[1:numsub,outcome],col=raincolors[dataframe[1:numsub,outcome]+1])

#}

1.8.2 The decorralted UMAP

  varlistcV <- names(varratio[varratio >= 0.025])
  topvars <- univariate_BinEnsemble(DEdataframe[,varlistcV],outcome)
  lso <- LASSO_MIN(formula(paste(outcome,"~.")),DEdataframe,family="binomial")
  topvars <- unique(c(names(topvars),lso$selectedfeatures))
  pander::pander(head(topvars))

X301_02, X520_38, X3325_02, X3439_22, X1752_66 and X4326_43

  varlistcV <- varlistcV[varlistcV != outcome]
  
#  DEdataframe[,outcome] <- as.numeric(DEdataframe[,outcome])
#if (nrow(dataframe) < 1000)
#{
  datasetframe.umap = umap(scale(DEdataframe[1:numsub,topvars]),n_components=2)
#  datasetframe.umap = umap(DEdataframe[1:numsub,varlistcV],n_components=2)
  plot(datasetframe.umap$layout,xlab="U1",ylab="U2",main="UMAP: After ILAA",t='n')
  text(datasetframe.umap$layout,labels=DEdataframe[1:numsub,outcome],col=raincolors[DEdataframe[1:numsub,outcome]+1])

#}

1.9 Univariate Analysis

1.9.1 Univariate

univarRAW <- uniRankVar(varlist,
               paste(outcome,"~1"),
               outcome,
               dataframe,
               rankingTest="AUC")

100 : X297_73 200 : X382_77 300 : X479_57 400 : X629_81 500 : X796_11
600 : X969_8 700 : X1173_74 800 : X1364_32 900 : X1599_26 1000 : X1848_42
1100 : X2201_1 1200 : X2498_63 1300 : X2860_34 1400 : X3257_2 1500 : X3669_6
1600 : X4208_08 1700 : X4836_61 1800 : X5379_65 1900 : X6237_78 2000 : X7115_25
2100 : X8025_19 2200 : X9106_13 2300 : X10504_79 2400 : X11887_21 2500 : X13358_42
2600 : X14760_46

univarDe <- uniRankVar(varlistc,
               paste(outcome,"~1"),
               outcome,
               DEdataframe,
               rankingTest="AUC",
               )

100 : X297_73 200 : X382_77 300 : X479_57 400 : X629_81 500 : X796_11
600 : La_X969_8 700 : La_X1173_74 800 : X1364_32 900 : X1599_26 1000 : X1848_42
1100 : X2201_1 1200 : X2498_63 1300 : X2860_34 1400 : X3257_2 1500 : X3669_6
1600 : X4208_08 1700 : X4836_61 1800 : X5379_65 1900 : La_X6237_78 2000 : X7115_25
2100 : X8025_19 2200 : X9106_13 2300 : X10504_79 2400 : X11887_21 2500 : X13358_42
2600 : X14760_46

1.9.2 Final Table

univariate_columns <- c("caseMean","caseStd","controlMean","controlStd","controlKSP","ROCAUC")

##top variables
topvar <- c(1:length(varlist)) <= TopVariables
tableRaw <- univarRAW$orderframe[topvar,univariate_columns]
pander::pander(tableRaw)

	caseMean	caseStd	controlMean	controlStd	controlKSP	ROCAUC
X301_02	1.0971	1.288	5.780	4.428	0.37901	0.901
X1714_44	0.9653	1.354	2.775	2.143	0.16574	0.876
X1752_66	0.4318	0.300	1.092	0.570	0.52591	0.868
X1539_08	0.4536	0.386	1.205	0.790	0.08336	0.866
X900_35	0.7958	0.696	2.402	2.230	0.03660	0.841
X440_87	1.8551	1.815	4.788	3.289	0.14184	0.839
X1522_71	0.5368	0.438	1.262	0.787	0.54259	0.836
X3325_02	39.9174	33.638	5.890	11.627	0.00130	0.833
X3428_53	0.0818	0.313	0.566	1.047	0.00325	0.829
X256_5	2.4015	2.773	8.994	9.942	0.00576	0.824

topLAvar <- univarDe$orderframe$Name[str_detect(univarDe$orderframe$Name,"La_")]
topLAvar <- unique(c(univarDe$orderframe$Name[topvar],topLAvar[1:as.integer(TopVariables/2)]))
finalTable <- univarDe$orderframe[topLAvar,univariate_columns]


pander::pander(finalTable)

	caseMean	caseStd	controlMean	controlStd	controlKSP	ROCAUC
X301_02	1.09708	1.288	5.780	4.428	0.379007	0.901
X1714_44	0.96530	1.354	2.775	2.143	0.165740	0.876
X1752_66	0.43183	0.300	1.092	0.570	0.525906	0.868
X1539_08	0.45356	0.386	1.205	0.790	0.083356	0.866
X900_35	0.79584	0.696	2.402	2.230	0.036597	0.841
X1522_71	0.53682	0.438	1.262	0.787	0.542593	0.836
X3325_02	39.91739	33.638	5.890	11.627	0.001299	0.833
X3428_53	0.08181	0.313	0.566	1.047	0.003253	0.829
X2870_03	0.65813	0.878	1.560	1.069	0.756294	0.821
X3439_22	24.04966	27.606	2.145	4.632	0.000556	0.816
La_X515_38	-0.00697	1.185	1.359	1.344	0.949004	0.805
La_X460_25	-0.08879	0.605	-0.779	1.147	0.309308	0.789
La_X440_87	0.25110	1.149	1.389	1.783	0.559404	0.747
La_X2671_39	0.51927	0.629	0.871	0.845	0.023880	0.722
La_X3540_88	-0.12276	0.775	0.103	0.198	0.810395	0.699

dc <- getLatentCoefficients(DEdataframe)
fscores <- attr(DEdataframe,"fscore")


pander::pander(c(mean=mean(sapply(dc,length)),total=length(dc),fraction=length(dc)/(ncol(dataframe)-1)))

mean	total	fraction
2.02	219	0.0825

theCharformulas <- attr(dc,"LatentCharFormulas")


finalTable <- rbind(finalTable,tableRaw[topvar[!(topvar %in% topLAvar)],univariate_columns])


orgnamez <- rownames(finalTable)
orgnamez <- str_remove_all(orgnamez,"La_")
finalTable$RAWAUC <- univarRAW$orderframe[orgnamez,"ROCAUC"]
finalTable$DecorFormula <- theCharformulas[rownames(finalTable)]
finalTable$fscores <- fscores[rownames(finalTable)]

Final_Columns <- c("DecorFormula","caseMean","caseStd","controlMean","controlStd","controlKSP","ROCAUC","RAWAUC","fscores")

finalTable <- finalTable[order(-finalTable$ROCAUC),]
pander::pander(finalTable[,Final_Columns])

	DecorFormula	caseMean	caseStd	controlMean	controlStd	controlKSP	ROCAUC	RAWAUC	fscores
X301_02	NA	1.09708	1.288	5.780	4.428	0.379007	0.901	0.901	0
X301_021	NA	1.09708	1.288	5.780	4.428	0.379007	0.901	NA	NA
X1714_44	NA	0.96530	1.354	2.775	2.143	0.165740	0.876	0.876	0
X1714_441	NA	0.96530	1.354	2.775	2.143	0.165740	0.876	NA	NA
X1752_66	NA	0.43183	0.300	1.092	0.570	0.525906	0.868	0.868	0
X1752_661	NA	0.43183	0.300	1.092	0.570	0.525906	0.868	NA	NA
X1539_08	NA	0.45356	0.386	1.205	0.790	0.083356	0.866	0.866	0
X1539_081	NA	0.45356	0.386	1.205	0.790	0.083356	0.866	NA	NA
X900_35	NA	0.79584	0.696	2.402	2.230	0.036597	0.841	0.841	0
X900_351	NA	0.79584	0.696	2.402	2.230	0.036597	0.841	NA	NA
X440_87	NA	1.85509	1.815	4.788	3.289	0.141837	0.839	0.839	NA
X1522_71	NA	0.53682	0.438	1.262	0.787	0.542593	0.836	0.836	0
X1522_711	NA	0.53682	0.438	1.262	0.787	0.542593	0.836	NA	NA
X3325_02	NA	39.91739	33.638	5.890	11.627	0.001299	0.833	0.833	5
X3325_021	NA	39.91739	33.638	5.890	11.627	0.001299	0.833	NA	NA
X3428_53	NA	0.08181	0.313	0.566	1.047	0.003253	0.829	0.829	0
X3428_531	NA	0.08181	0.313	0.566	1.047	0.003253	0.829	NA	NA
X256_5	NA	2.40151	2.773	8.994	9.942	0.005764	0.824	0.824	NA
X2870_03	NA	0.65813	0.878	1.560	1.069	0.756294	0.821	0.821	0
X3439_22	NA	24.04966	27.606	2.145	4.632	0.000556	0.816	0.816	1
La_X515_38	- (0.157)X308_56 + X515_38	-0.00697	1.185	1.359	1.344	0.949004	0.805	0.814	-1
La_X460_25	+ X460_25 - (1.208)X681_12	-0.08879	0.605	-0.779	1.147	0.309308	0.789	0.371	-1
La_X440_87	+ X440_87 - (0.177)X629_81	0.25110	1.149	1.389	1.783	0.559404	0.747	0.839	-1
La_X2671_39	- (6.958)X877_66 + X2671_39	0.51927	0.629	0.871	0.845	0.023880	0.722	0.699	-1
La_X3540_88	- (0.044)X3325_02 + X3540_88	-0.12276	0.775	0.103	0.198	0.810395	0.699	0.758	-1

1.10 Comparing ILAA vs PCA vs EFA

1.10.1 PCA

featuresnames <- colnames(dataframe)[colnames(dataframe) != outcome]
pc <- prcomp(dataframe[,iscontinous],center = TRUE,scale. = TRUE)   #principal components
predPCA <- predict(pc,dataframe[,iscontinous])
PCAdataframe <- as.data.frame(cbind(predPCA,dataframe[,!iscontinous]))
colnames(PCAdataframe) <- c(colnames(predPCA),colnames(dataframe)[!iscontinous]) 
#plot(PCAdataframe[,colnames(PCAdataframe)!=outcome],col=dataframe[,outcome],cex=0.65,cex.lab=0.5,cex.axis=0.75,cex.sub=0.5,cex.main=0.75)

#pander::pander(pc$rotation)


PCACor <- cor(PCAdataframe[,colnames(PCAdataframe) != outcome])


  gplots::heatmap.2(abs(PCACor),
                    trace = "none",
  #                  scale = "row",
                    mar = c(5,5),
                    col=rev(heat.colors(5)),
                    main = "PCA Correlation",
                    cexRow = 0.5,
                    cexCol = 0.5,
                     srtCol=45,
                     srtRow= -45,
                    key.title=NA,
                    key.xlab="Pearson Correlation",
                    xlab="Feature", ylab="Feature")

1.10.2 EFA

EFAdataframe <- dataframeScaled

if (length(iscontinous) < 2000)
{
  topred <- min(length(iscontinous),nrow(dataframeScaled),ncol(predPCA)/2)
  if (topred < 2) topred <- 2
  
  uls <- fa(dataframeScaled[,iscontinous],nfactors=topred,rotate="varimax",warnings=FALSE)  # EFA analysis
  predEFA <- predict(uls,dataframeScaled[,iscontinous])
  EFAdataframe <- as.data.frame(cbind(predEFA,dataframeScaled[,!iscontinous]))
  colnames(EFAdataframe) <- c(colnames(predEFA),colnames(dataframeScaled)[!iscontinous]) 


  
  EFACor <- cor(EFAdataframe[,colnames(EFAdataframe) != outcome])
  
  
    gplots::heatmap.2(abs(EFACor),
                      trace = "none",
    #                  scale = "row",
                      mar = c(5,5),
                      col=rev(heat.colors(5)),
                      main = "EFA Correlation",
                      cexRow = 0.5,
                      cexCol = 0.5,
                       srtCol=45,
                       srtRow= -45,
                      key.title=NA,
                      key.xlab="Pearson Correlation",
                      xlab="Feature", ylab="Feature")
}

1.11 Effect on CAR modeling

par(op)
par(xpd = TRUE)
dataframe[,outcome] <- factor(dataframe[,outcome])
rawmodel <- rpart(paste(outcome,"~."),dataframe,control=rpart.control(maxdepth=3))
pr <- predict(rawmodel,dataframe,type = "class")

  ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
  if (length(unique(pr))>1)
  {
    plot(rawmodel,main="Raw",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
    text(rawmodel, use.n = TRUE,cex=0.75)
    ptab <- epiR::epi.tests(table(pr==0,dataframe[,outcome]==0))
  }

pander::pander(table(dataframe[,outcome],pr))

	0	1
0	31	6
1	4	105

pander::pander(ptab$detail[c(5,3,4,6),])

	statistic	est	lower	upper
5	diag.ac	0.932	0.878	0.967
3	se	0.963	0.909	0.990
4	sp	0.838	0.680	0.938
6	diag.or	135.625	35.975	511.303

par(op)
par(xpd = TRUE)
DEdataframe[,outcome] <- factor(DEdataframe[,outcome])
IDeAmodel <- rpart(paste(outcome,"~."),DEdataframe[,c(outcome,varlistcV)],control=rpart.control(maxdepth=3))
pr <- predict(IDeAmodel,DEdataframe,type = "class")

  ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
  if (length(unique(pr))>1)
  {
    plot(IDeAmodel,main="ILAA",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
    text(IDeAmodel, use.n = TRUE,cex=0.75)
    ptab <- epiR::epi.tests(table(pr==0,DEdataframe[,outcome]==0))
  }

pander::pander(table(DEdataframe[,outcome],pr))

	0	1
0	31	6
1	3	106

pander::pander(ptab$detail[c(5,3,4,6),])

	statistic	est	lower	upper
5	diag.ac	0.938	0.886	0.971
3	se	0.972	0.922	0.994
4	sp	0.838	0.680	0.938
6	diag.or	182.556	43.143	772.462

par(op)
par(xpd = TRUE)
PCAdataframe[,outcome] <- factor(PCAdataframe[,outcome])
PCAmodel <- rpart(paste(outcome,"~."),PCAdataframe,control=rpart.control(maxdepth=3))
pr <- predict(PCAmodel,PCAdataframe,type = "class")
ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
if (length(unique(pr))>1)
{
  plot(PCAmodel,main="PCA",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
  text(PCAmodel, use.n = TRUE,cex=0.75)
  ptab <- epiR::epi.tests(table(pr==0,PCAdataframe[,outcome]==0))
}

pander::pander(table(PCAdataframe[,outcome],pr))

	0	1
0	33	4
1	11	98

pander::pander(ptab$detail[c(5,3,4,6),])

	statistic	est	lower	upper
5	diag.ac	0.897	0.836	0.941
3	se	0.899	0.827	0.949
4	sp	0.892	0.746	0.970
6	diag.or	73.500	21.908	246.593

par(op)

1.11.1 EFA

  EFAdataframe[,outcome] <- factor(EFAdataframe[,outcome])
  EFAmodel <- rpart(paste(outcome,"~."),EFAdataframe,control=rpart.control(maxdepth=3))
  pr <- predict(EFAmodel,EFAdataframe,type = "class")
  
  ptab <- list(er="Error",detail=matrix(nrow=6,ncol=1))
  if (length(unique(pr))>1)
  {
    plot(EFAmodel,main="EFA",branch=0.5,uniform = TRUE,compress = TRUE,margin=0.1)
    text(EFAmodel, use.n = TRUE,cex=0.75)
    ptab <- epiR::epi.tests(table(pr==0,EFAdataframe[,outcome]==0))
  }

  pander::pander(table(EFAdataframe[,outcome],pr))

	0	1
0	31	6
1	4	105

  pander::pander(ptab$detail[c(5,3,4,6),])

	statistic	est	lower	upper
5	diag.ac	0.932	0.878	0.967
3	se	0.963	0.909	0.990
4	sp	0.838	0.680	0.938
6	diag.or	135.625	35.975	511.303

  par(op)